{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T03:28:20.690571Z",
     "start_time": "2019-07-13T03:28:20.687508Z"
    }
   },
   "source": [
    "This notebook demonstrates how we can build a working linear regression model using NumPy library only\n",
    "\n",
    "Assume we know X and y are linearly correlated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Initialize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T05:53:04.001272Z",
     "start_time": "2019-07-13T05:53:02.607674Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "np.random.seed(42) # for reproducibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Generate data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure X and y are linearly correlated with normally distributed noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T05:53:30.658718Z",
     "start_time": "2019-07-13T05:53:30.655157Z"
    }
   },
   "outputs": [],
   "source": [
    "n = 10000 # 10000 datum points\n",
    "X = np.random.uniform(-10,10, n) \n",
    "noise = np.random.normal(0, 3, n) # Gaussian distribution\n",
    "true_w, true_b = 7.6, -3.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T05:53:59.353611Z",
     "start_time": "2019-07-13T05:53:59.350757Z"
    }
   },
   "outputs": [],
   "source": [
    "y = true_w * X + true_b + noise # y = w * x + b + ε"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T05:54:06.347453Z",
     "start_time": "2019-07-13T05:54:06.134055Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xlc1HX+wPHXZxiGS+QSUbxANI1KTfHItDTttqw27c5OsrZz292Odbdzd9vOn9VW2qltW2r3urWlZqaZJpqZkgQjqHjgwCDIcAwDn98fM/N1gBnAOER4Px8PH/D5zne+ny/D+OYz7+/n+/4orTVCCCE6P9PRPgEhhBDtQwK+EEJ0ERLwhRCii5CAL4QQXYQEfCGE6CIk4AshRBchAV8IIboICfhCCNFFSMAXQoguwny0T8BXjx49dFJS0tE+DSGEOKZs3LixUGsd39R+HSrgJyUlkZGRcbRPQwghjilKqZ3N2U9SOkII0UVIwBdCiC5CAr4QQnQREvCFEKKLkIAvhBBdhAR8IYToIiTgCyFEFyEBXwghjiK7w8m8VVbsDmeb9yUBXwghjhK7w8m9izfz98+3c+/izW0e9CXgCyHEUbIkYzcrs2ykxEewMsvGkozdbdpfi0srKKWGAIt8Ng0E/gJEAzcDNs/2B7XWn7W0PyGE6CxGJ8WSEh/B3VMG8+EPe5iamtCm/bU44Guts4ARAEqpIGAP8BFwPfCc1vrplvYhhBCdhd3hZEnGbkYnxXLDWxs4WFHNnE+2UlLhYljfvdxz5nFt1ndrF0+bAli11juVUq18aCGEOPZ4A/yMtH7ERlhYsDaXuSty6BsdxsGKagBKKlwAVDhdbXourR3wLwfe9WnfrpS6FsgA7tVaF9d/glIqHUgH6N+/fyufjhBCHF1LMnbz98+3U1RWxY/5B8m1OQCwHarws3fbDpRb7aKtUsoCXAgs8Wx6GUjBne7ZBzzj73la6/la6zStdVp8fJPlnIUQ4pjgnW45OimWkf2jWPDdTtbnFnOgzD0Tp6rGz5PaODHSmiP8c4FNWusCAO9XAKXUq8DSVuxLCCE6FG/qZmpqAsszC9hTXMHCdTsJCVJU1eiAz4uscvDIsld4ffRFoAe26Tm2ZsC/Ap90jlKqt9Z6n6d5MbC1FfsSQoh2Vz8f72W1lZG+MAOrzcE3v9j41lpEv5gwgMDBXmue+N8LXL7lSwByYxKpYmqbnn+rpHSUUuHAmcCHPpufVEr9pJTaAkwG7mmNvoQQ4mjx5uPrz5d/6JOtWG0O+sWE4ap1B/hTUuIwB0jRnJm9jrwnLzCC/Stjf8MLp15BW+d0WmWEr7UuB+LqbbumNY4thBAdxYy0fnW+gnt0v8NzIbagtJLdxRX0iwlj+74SXPUG9wmHCln/0nVGO797T6be9BKVwaEAVFb7S+y3ng61pq0QQnRksREWZqT1q5OrX51dyN6SSgCcnvTN7uIKdvvMSTTV1vCvRX9m/K4txrazbniRX+KT6hw/q6C0Tc9fAr4QQhwBb1pn3Y4iVmbZCA5qPA1z9Q+f8fiXLxntOWfeyr9Gnu933+F9Y1r1XOuTgC+EEEdgRlo/yp0uKqpr+c5aSGX9vI3HcbY8vnzjdqO9tv8wrr7sMWpNQXX2U4AGosLMzJ6U0oZnLgFfCCGaZLWV8fjSTOZMSwXgw0172F1cQXJcOLlF5XX2DamuYsVrt9K39ICxbdytb7G/ew+/x/b+uZg+vE+dmT9tQQK+EEL4YXc4WbA2F1Cs21HI+txiSip+JLfQQXG5uyRC/WB/39dvcev69412+sV/4svjTmmyr34xYdzdhjV0vCTgCyGEHwvW5jF3RQ4AwZ4J7D/vK6WiurbBvuN2beG9dx802u8NO4v7z7kDmllT7NwTe7X56B4k4AshhF/FjioALEHKmH1TP9jHlJfwwwtXGW1HcCjjfruAQyERjR67V/cQekeFclKfaGIigpk1PrmVz94/CfhCCFGP3eFkZZZ7KQ+z6XDAN2jNi58+ybTtq41Nl1z9FJv6HB/wmN6LswAD4iJ4+epR7TKq9yUBXwjRpfnm6meNTwLg1n9tZHexu5pleb1R/YWZq3j+P08Z7acnXs2L4y9vsp9rxvXnm+xC8orKWZ9rZ0nGbm45vW1n5dQnAV8I0WXZHU7ufHcTa3KKANi4s5jU3t1Zn2tvsG+/g/tZPe8mo/1zfBLTr30Opzm4zn6+I/mU+AjsDidnHt+Tu88cwt1nDmHB2jxA17lbt71IwBdCdFmvfJ1jBHuANTmFZO49WGcfc42LD/71B4bvzza2Tbp5Hnmxffwe0zf5k28vp6pG832e+7bb2AhLm65o1RQJ+EKILsPucPLK1zn8mF/CkIRuRp7epMBT8wx7+eFVp9LXf8CDX79ptO897x4+OGlKs/qKCDHhqHKng/KKyo9KCqc+CfhCiC5jScZu5q/OBTDSNiFmE1Wuunn6k/Zl85+Fhwv8fjl4HLdc/CBa+S8wHGxSVHv+YphNClet5sTEaMYNjKPCWUOYxXRUUjj1ScAXQnQ6gerWT01N4KWvcyipcGEGXFAn2EdUlbPmlRuJqTxkbEu7/W0KI/zXuAkOgvhuofzt4pP4KusAq7Js/P6s4/jwhz3MmZZKSny3Rs+nvUnAF0J0Gt7AWu6sYe6KbP5v+S88+Zth5NjKqKiuJSPPTkmFi+AgRXW9qZaPLHuZWZv+a7Svmfkoq5NHNtpfdQ3sLanktTW5TBzcg532cvaWVPLm9WPq7OctuAYc1bSOBHwhRKfgu+rUXVMGERpsoqK6lt+/v6VBysY32E+ybuCt9x8x2q+nTeexKTcH7CdIQY1235AVbDLhqK4htXek31r5Xo091p4k4Ashjlm+qZLHl2ZitTlIiY9g1vhk9h2sZPHGfFLiI8jcd6jBc+PLitnwz8PrNB2IiGFS+nzKLWEB+zN5gj24a9/fOmkg4ZYgI1UTaPTe2GPtSQK+EOKY5ZsqueOMweQWOhg1IIZXvs4hOtxC+mkDySk4VCfgK13LG+8/wuQdG41t5183l20JgQOydyFy70yexKhQpg1PZNb4pKOakz9SEvCFEMcsb4pkdFIsd733A7uLK8jzqWDZs5uFA2VOoz3zxy958n/PG+1Hz7iZN0ZPb7KflJ7dmDAoHoDMfSWsySkiLsJyTAV7kIAvhDhGuUsi5FHhdHHPos1GKQQAswlctRjBPqVoNyteu9V4PKPP8Vx25RPU1FuMpD7vyD5z3yGmj+jDLaen1EkjHWtaLeArpfKAQ0AN4NJapymlYoFFQBKQB8zUWhcHOoYQQvhTf1qj3eHk3sWbjRun6vNeow1xOfnszTtJsecbj506+w32RPVsVr9VNZqkuHCmj+hjBPiOko//NVp7hD9Za13o074fWKG1fkIpdb+nfV8r9ymE6OS8ufpvfrGRlhRLsaeaZbBJYQ7SVNccDvImoBa4Z/U73LX2XeMYt06/n8+HTmh2nyP7RxNuCeKR6Sca8+mPdW2d0pkOTPJ8vwD4Ggn4QogmeEf0o5NieeGrbO44YzCTh8SzMsvGt9YiIkLcd7xW12rqr0dycn4mH7zzR6P94QmT+d35v2t0MRLl+VeL+4LsjLR+RuVMb7GzWeOTj7mcfX2tGfA18KVSSgPztNbzgQSt9T4ArfU+pVTzPkcJIbo074g+NiIYu6OaXfZynrp0OCUV1fy8t8SoUeOre2UZGS9cjaXWXQvHpUyMuuMdSsIiA/bTLSSIEf2iWZNTZBQ9mzY80ShwNm+Vlbkr3EXTwi3mYzaV49WaAf9UrfVeT1BfppTa3pwnKaXSgXSA/v37t+LpCCGOFd4R/dTUBJZnFjA1NYHV2YWsySnEpMBqczBz3ne4anXDJ2vNM/99lt9sW2lsmnHlE2zod2KT/Z53Yi9umTSIG97cwE57ORMG9WC2T1CfkdaPcmcNR6uccWtrtYCvtd7r+XpAKfURMAYoUEr19ozuewMH/DxvPjAfIC0tzc9vUwjR2XlH9PO+sWJ3uBcIT+0dyZqcQmPuu79gf07Wt7zy8d+N9tzxl/PcxKub3W90RAjLMwvYaS9n8pB4npk5ok7a5miXM25trRLwlVIRgElrfcjz/VnAo8CnwCzgCc/XT1qjPyFE5zIjrR+LM3ZjtTnoGx3GF9v2k28vD7h/YukB1r58g9HeEZPIuTe8SJW56Rx7uMVEdJiFvSWVhAWb6pQ9ONZz9E1prRF+AvCRcl8UMQP/1lr/Tym1AVislLoR2AXMaKX+hBCdhHc+/YmJkdgdTsJDTGzaddDvvkG1Nbz37wcYvSfT2Dblxpex9vCfbvHO2Ak2KWIjLBQcqqLcWUu5s5LJQ+K5cESfDlHFsr20SsDXWu8AhvvZXgQ0b7UAIUSXsWlnMfcs2szYgbHsPVhRZ9Wp0opqv8+5LuNTHl4x32jfd84dLBp+doP9FDAwPgKnq9a4GWtwzwgiwyz0iQll064SkuLCmTMtleWZBR2iimV7kTtthRDt7g/v/8hOezk7PWmbIJ+iZPWqFnP8gR18/uadRntV8kium/FwwMVINO6LvABJceEkdA9hfa77fs+xyTHG9M7lmQUdpople5GAL4RoV1ZbGVFhwfTuHkp8pIUte0obBHmAMGclX7+aTkLZ4QXFx9y2gAORcU324b5pysyanEKmj0gEFOtz7QzvG8PsSSl10jhdYWTvJQFfCNGmrLYyHl+ayZxpqcSEW7jxrQ3kFZWTFBfO784cwp8/2cqhymoOVhxeS3bOile5KePwHI/rLn2Ir1NGN6u/pLhwXpvl3tcb2GeNT+6yQd6XBHwhRKsItIzfQ59sZU1OEdU1Wxk1INaoZplXVM79H25hf2mVse+peZt5Z9Eco/32yefx57Nua7TfXt1DCDIpokKDiAwL4W+XnGT07xvYu2qQ9yUBXwjRKrxz6cudLsItZiPwJ0a7FxTZtNPO9n2ldZ7jDfZxjoNsfPHw/PmDod2YMPsNykLCm+x3QFwE63PtXDplEOEWMzHhnX+2za8lAV8I8av5juq9Fz7LnTWewF9DhbOGz37a595erSmvrjsDR+la5n30N87KXmdsu+Da5/ip9+Am++4XE0ZidBj3nTOUDXl2o1+Q0XwgEvCFEL9a/cW5Z6T1Y8HaXGaO6stb3+ZSUuny+7wgBdN/WsGz/33O2Pb3Sdcxb+ylzeo3MSqUS0b2Ze6KbDbk2Y069d7lBoV/EvCFEL+at+ZNUVkVVlsZD32yjTU5hcSEBwcM9kn2PXz96i1Ge0uvQVxy9dO4gvyHo2CTYnBCN4od1ewrrSQpLpzXrxtNTLilToDvyhdjm0sCvhCi2eoXOSsqc7Imx13kbMX2A8b890OV1YSaFZWuw/Mtg2uq+XTBPRxvyzO2TbzlNXZH92q0z+pabaxJe2pKHC9cORKgS90h21ok4Ashmm3B2jzmrshmdbaNNTlFTBjknhMfFWrGanMYSwK6ausWO7vtu8X88ZuFRvvOC/7Ap6mnB+wnOiyIgxU1DbanJcUSG2Fh3iqr5Ot/BQn4QogmeefS94txz7hJ7R1Fau8o1mTbCA5SRvqmqt4dVCP2ZvHx2/ca7aVDJ3L7hX9sdDGSfjFhlNQrrzAgNpyzT+hlLErS1e6QbS0S8IUQAXlTOF9tP8D6XDtjk2N54NyhTE1NIH1hhpHCqS+yysG6f84iorrS2HbyHe9QHB4VsK+QIBNx3SzsLq7gnBMS+DKzgFrt/gOw015OXDeLkb6RfP2vIwFfCBGQdxbO2OQYwH2T07vf7+KLbfux2hwooM6YXmv+/r8XuGLLl8amKy7/K98NaFBbsQ53KqiWvSXuKpZlVTVGHfxzT+xNXDeLjOZbgQR8IUQdvhdmy501pE9MBhQD4iL4YGM+Ndp9l2z9YD81ez2vffiY0X5l7G94YtL1zerzstH9iYmwAJoLR/Rh0fe7qK6pZXjfKGZPSiE2woLd4WTeKqtcqG0BCfhCiDo3UHkvzH6xbT+bdh2ke6iZ0koX3UPNdYqceb9NOFTI+peuM7bnd+/J1JteojI4tNE+lQKtYfKQeO4+8zgjiM9bZWX+6lweOHdonbRN/Tn/4shJwBdCGMF0z8EKPtm8B8BYhKTUc0HWWV13Xr2ptoa3F/+ZU3duMbadfcOLZMUnNdpXZIiJQ1W1XDisF6mJ0Q1G7IEuyMqF2pZTWnecZWTT0tJ0RkbG0T4NIboM7+ybO84YzKpfbI3eHevrqh8+469fvmS055x5K/8aeX6jz0ntFUlkmJkhvboTEx7MrPHJkpppJUqpjVrrtKb2kxG+EF2U3eE0ZtrkFjroEx3WZLA/zpbHl2/cbrTX9h/G1Zc9Rq0pqNHnjU2OZdzAOOauyGZ9bjGTh8S3ys8gjowEfCG6GKutjAc//Im9ByvYXVxBqFmRV1ROXlE5kSFBHKpqeMNTSHUVy1+/jX4lBca2U259k33dmxe4zSbF6cfFs3GnneqaWlZm2ViSsVty8e1MAr4QXczjSzNZn+teRSoqzEyJz8IjmtoG+/9h1QJ+u26J0U6/+E98edwpjfbhu2RhUlw431qLsHyVzZqcIu6aMpgzhiZILv4okIAvRBcz65QkMnbaCTGbKK+sO5ovqzp8TW/cri289+6DRvu9YWdx/zl3NHqXLMBxCRFMOi4B0IRZzFw4IpHlmQVMTU1g3MACmVZ5FLU44Cul+gELgV5ALTBfaz1XKfUwcDNg8+z6oNb6s5b2J4RoHn8rUNkdTs+SgjUcomHqBiC6opTNz19ptMuDQxh32wJKQ7s1q1+nSzN/9Q4mD4nnmZkjiI2wkHK6+7ner+LoaI0Rvgu4V2u9SSkVCWxUSi3zPPac1vrpVuhDCNEE37VjU+K71Zm3PjU1gQc//IlcWxkHypz+D6A1L3z6JBdsX21suuTqp9jU5/hG+w0xK6o8VTH7RIfy7MwRvPBVtuTpO6AWB3yt9T5gn+f7Q0qpn4E+LT2uEKJpvqP4x5dmsjLLBmTyzMwRlDtruGvKoCbr3gBckLmKF/7zlNF+euLVvDj+8madw/Thfdhpd7A+t5jzT0pk5IAYnpk5wjgv0XG0ag5fKZUEnAysB04FbldKXQtk4P4UUOznOelAOkD//v1b83SE6PR8R/FzpqVS7txCSUU1t/4rg/W5xYQFm1iWWRAw2Pc9uJ81824y2j/HJzH92udwmoObfQ69o8PoHR3m7s/inp4pxc06plYL+EqpbsAHwN1a61Kl1MvAY7jvwH4MeAa4of7ztNbzgfngvvGqtc5HiK5gRlo/yp0uyp3ufHxhmdMI7iYFFdW1xuIhvsw1Lt5/5w+M2JdtbJt88zxyY5v/4bxfTBiXjOxrlCyW5QU7vlYJ+EqpYNzB/h2t9YcAWusCn8dfBZa2Rl9CiMNiIyyEW8z8/fPtzP/GSkV1LSFmE1WuWmo1KCDIpOosRnLz+g/509dvGO17z7uHD06a4vf4wSao9pmpmRQXzvC+0fyYf5BnZ45g5IAY4zEZ0Xd8rTFLRwGvAz9rrZ/12d7bk98HuBjY2tK+hBCHWW1l/OXjraT07Ea3kCDKqmpQgNkEVZ59NIdXnjpxfw5LF9xtPP/LweO45eIH0crk9/i9uodw1gm9WPjdTsAd7POKyknuEUFeUTkb8ux1Ar7o+FpjhH8qcA3wk1Jqs2fbg8AVSqkRuN9zecAt/p8uhDhSdoeTG97cwE57Od9ai/DOjNeAw1n35qmIqnJWz7uJ2IpSY1va7W9TGNF4sD4rtRehZhNjk2MZ3jeKy8b0bzCfXhxbWmOWzhrA350YMudeiFZkdzhZsDaPCqeLH/NL2GkvBxqmXXzvcn1k2cvM2vRf47FrZj7K6uSRDY5tUpDcI4LdRQ6ctTCyfzQxEcHMXZEDwBlDe5IS303m0x/j5E5bIToAfzdJ1bckYzdzV2T7eaTuUiQ1Gk7fsZEFSx4ytr2eNp3HptwcsP9aTZ2ZPBMH92DW+GTj+DKa7xwk4AvRAdRf3MN31alPN++h2FHN1r0lDOvTHVdNLSWVLpzVtdgcTqp9LsjGlxWz4Z/XGG1bRDSnp79KuSUsYN/hwSbKq2uJCFFEhoZw2qA4vB/a7zlzSNv8wOKokIAvRAcwNTWBdTuKmJqaABz+A7A6u5A1OYV+n2PySaQqXcsb7z/C5B0bjW3nXzeXbQmBZ854PxdcfUoSK352z9V3VFWyt6SSxRv3EG4Jkpk3nYz/y/NCiHa1PLOAlVk2lmcWYHc4KXfWMHNUX37ac7DOfmaf/7Hegf2MLV+S++SFRrB/9IybSbpvaaPB3oQ72I9JiiEsOIinLh1OUlw4AKm9o3jg3KGSxumEZIQvRAfgDa6jk2K55KVvySsqN9aS9Z2B4/K5OJtStJsVr91qtDcmDmXmVf+gppHFSCYOisMcZKJfTDgL1+3kwKEq5q7IJtwSxIe3ndrkdQRxbJOAL8RRUP8irbcUwfVvfk9ekXv2jXct2fq3n4e4nHz25p2k2PONbafOfoM9UT0D9uedyZPcoxuPXnQiVlsZ31oLsdocTB4Sb5yHpHA6Nwn4QhwF9S/Sfr39AHe+9wPBpsZrzd+15t/c8+2/jfat0+/n86ETmuzPO21zR2EZ4E4heYO9t4Sx6Pwk4AtxFHhTOFNTE/jLxz/x9rpdxkg+1KyodNUd14/Kz+SDd/5otD88YTK/O/93TS5G4nVcQgQ9I0N5ZPqJdfqX9E3XIgFfiDbmO8XSe6fq8kz3napLMnazcN2uOvv7BvvulWVkvHA1llp3eselTIy64x1KwiID9hdhMdW52/bUlDheuHJkncAu6ZuuSQK+EG3Mm75Zt6OIlVk25n1jxe6optzpothRTYQlCIez3upTWvPMZ8/xm61fGZtmXPkEG/qd2GR/DmetsSjJ2OSYBsEemnejl+h8JOAL0cbcJYxrqHC6yDlQxu7iCiJDg3hn3S4KHQ1Xnzo7ay3zPv6b0X7+lMt49rRrGuznKyQIqjx/MwbEhrPTXt5ofr7+NQTRNUjAF6IV+Rs55xU6eHtdHnZHNb26hwD4XVM2sfQAa18+vGTEjphEzr3hRarMTY/Aq2ogMSqUgfHduGlCMq+u3sHghMBpH98cvug6JOAL0QL1A7x35FzkcJJdcIhZpyQx+18bqXTVEqRgf2lVg2ME1dbw7rsPMCY/09g25caXsfZoOhhHhZkpqXDn9/vFhrEmp5DgIMW31iK+tRYRFyBXLzn8rkkCvhAtUD814h0xf7W9gPW5xazOLjTq0df4Wc9t1sb/8MjyeUb7/rNv570R5wTsz2yCHt1CiA230C3UzG2TBvHq6h2c0CeKy0b3My4KD+u7Byl6JuqTgC9EC/hLjRQ5nOTa3DdPuWp1vVqWbscf2MHnb95ptFclj+S6GQ8HXIzEKzQ4iP2lVewvreKBc4eSVXCIb61FnHZcfJ3yxVL0TPgjAV+II2S1lfH40kzmTEslJb6bMb1yamoCD374E+tz7XX29w32Yc5Kvn41nYSyw/uMuW0BByLjmtV3WVUNp6bEkZYUW+ePjIzkRXNIwBfiCD30yVbW5BRRXbOVf900jgVrc5m7Iof/btnLlj2lAZ/3p69e4+YNHxvt6y99iJUpoxvs5+8TgVdUmLnBNEvJxYvmkmqZQjTB7nAyb5UVu2cKZWrvqDpfKzw3OQUK9qfmbSbvH9OMYP/2yeeRdN9Sv8Ee4OwT3CWSu4ceHo+FW9x31E4f3kfmzYtfTUb4QjTB98Ls1NQEMveVkD5xIONT4pj01EosZv/lDeIcB9n44tVG+2BoNybMfoOykPBG+/P+YRnUM4IT+0RjPVDGvWcNYUOeXVI3okUk4AtRT/2plr4XZu9dvNmTzqnl7XV5VFTXNjyA1sz76K+cnb3O2HThtc+ypfdxTfbdLyaMob27831eMZt2lRAVZnFflM2zS+pGtFibB3yl1DnAXCAIeE1r/URb9ynEr+EN9OXOGuauyKbcKHeguXBEH5Zk7OYsz8pUhYeq/Ab7i7d+xXP/fdZo/33Sdcwbe2mj/YYFm3DWaAbFR5BVUMYOm4P0iQMJs5i4cEQfxg0skJG9aBVtGvCVUkHAP4EzgXxgg1LqU611ZuPPFKL9eVM36RMHMnlIPMUOJwvX7QRg2bZ9ZO53EKTc8+mtheV1njugeC+r5qcb7S29BnHJ1U/jCmr8v5jZpLhgWG8Wb9xDj24hRIdbWJNTyKgBMdxzpvsTgXeqpRAt1dYj/DFAjtZ6B4BS6j1gOiABXxw1gQqHeUfR5U4XK7NsxpJ/AD/vdwANb54Krqnm0wX3cLwtz9g28ZbX2B3dK2D/KT3CjT8YrlpN7+hwJg+JZ2WWjVNTvNMzA83TEeLXa+uA3wfY7dPOB8a2cZ9CNKr+3bG+5YsBTj+uJ0u37MNqc5AYFcrekkq/4fe27xbzx28WGu07L/g9n6ZOarTvqDAzPSJDsBa6lzC8fEx/Zo1PMs7Lt3SyEK2trQO+v+kLdf7vKKXSgXSA/v37t/HpCNHw7ljvH4DV2YWsySlkwqA4rDYHJqCgtBKouyjJiL1ZfPz2vcbxlg6ZwO3T72vWYiTTh/dh1qlJpC/MwGpzEBYcZHza8F6UlRSOaCttHfDzAd+hSl9gr+8OWuv5wHyAtLQ0+Rwr2lz9wmEz0vpRVObks5/2AVBaUQ1ALRjDk0qXJrLKwXcvXUc3Z4Xx3JPveIfi8Cij3TPCgs3hbPCJYGT/KCYO7sms8UnuImuzx3suELukTLFoN20d8DcAg5VSycAe4HLgyjbuU4gm1V+FKnNfCfkHKwgJUg1voNKav33xIlf++IWx6YrL/8p3A4Y3OG5hecNgHxZs4rVZY4zrBb7XEADCLWZJ4Yh20aYBX2vtUkrdDnyBe1rmG1rrbW3ZpxD+eIPs6KRYXvgqm/huISzemM/cFVmUOzVJsWEEKaiqd1V2Ss56Xv/gMaP9ytjf8MSk6wP2U6vdefpQcxAFh6owmxQvXzWK2AiLz7RPF3NX5ADuUb2M7EV7afN5+Frrz4DP2rofIRrjrXfTPdRMaaWLIE9RkXKnO8Dn2Svq7J9wqJD1L11ntPO792TqTS9RGRzq9/gmPCkg3Hn6i07uwx/e/5H/58fMAAAgAElEQVSnLh3OyAExwOFrBXdNGcwD5w6VUb1od3Knreh0/E279Na7Ka10LxZSU4sxp96XqbaGhYv/woSdPxrbzr7hRbLikwL2ZwlSOGs0MeHBFJdXExNhYUOeHavNwYY8uxHwfS8WSz0ccTRIwBedjr9pl+t2FAJ1g3z9YH/l5s/52xf/NNp/PnM2b4+c1mR/zhpNbEQwr107mmeX/UKF08WFIxKBumWLZZUpcbRJwBedhm+efsKgOIrKnFhtZTz0yTbjQqy/VacG23ay7I3fGu21/Ydx9WWPUWsKala/MeHuYL8hz86aHPfUzrhuIRLcRYcjAV90Gt6R/eQh8azJKWJNThGLMnYZa77WF1JdxbLXb6N/SYGx7ZRb32Rf9/hG+wkOUlR7/nJMGBTH81e469Mn9Yjw1N/Rkp8XHZIEfHFM883XT/UUNpt1ShI7bA522ssDBvvff7OQ279bbLTTL/4TXx53SrP6PL5XJFv2lHJqyuFgD+6Ujbf+jRAdkQR8ccyyO5zcu3gzK7NsfPOLjYrqGjbtOsjP+0rZX1rl9zljd/3EoncfMNrvDTuL+8+5I+BdshGWIHpFhWC1HS6WNi6lB+cPS2RqaoLfmjxCdFQS8MUxye5wcue7m1iTU0T3UDPfWouMx/wF++iKUjY/f/ievwpzCGN/u4DS0MBlDIKDFG/fOJakHhE88dnPfJG5n7NTezH79BRiIyzMW2WVu2TFMUUCvujw7A4nC9bmAZpZ45MBuPPdH1iT4w7y3qmWfmnN8/95igt//sbYdMnVT7Gpz/FN9ltdo9mQZ2dDnp3FG/MBSOnZzRjN16/JI0RHJwFfdCj159D7pm3AXYYAYE1OYZPHmvbzN7z46ZNG++mJV/Pi+MsbfU5YsMlY2GRscmydksmgZJqlOKZJwBcdiu8ceu+Sgt7a9KcNjqeorIpKVy2RIUEcqqrxe4y+B/ezZt5NRjurR38umDUXpzm4zn6KhlXno8KCqaiuYmxyDC9fPcoYzd9z5pBW+xmFOFok4IsOxTdNsiRjNyuzbESFmckrKqem9gC7iysCPtdc4+L9d/7IiH2/GNsm3zyP3Ng+fvfXQM9uFg6UOZk+vDellS7uOGOwsVi4XIgVnY0EfNFh1K8iWe50MWFQnJGrbyzY3/T9h8xZ+YbRvve8e/jgpClN9lnpcqdvihzV/Osm99o83lIIQnQ2EvBFh+GbzvFWlEyfOJAdNgd7Syr9PufE/TksXXC30V42aAzpl8xBK1Od/bzFzbxpnDCzosKlGdanO/tKq/idzJ8XXYAEfNHu/BU3szucFJU5GZscQ5HDaeybua+ECYPiWLxxT51jhDsrWP3KjcRVHK5dn3b72xRG+B+deytZXuhJ3XgvAm/bd4ji8mo25NlJ6hEh8+pFpyYBX7Sbw/Xga5i7Ihtw5+oXrM1l486Dxsyb9bnF9OoeQoTFxJqcIkz1jvPQ8nlcv/E/RvuamY+yOnlks85hf2klw/tGEx8ZwnfWIk4ZGEvv6HDjmoHMqxedmQR80W68AXVMUgwTBvVgamqCUaceoGdkCAcOuW+a8r15yjs6P33HRhYsecjY/uaoC3hk6i0B+1NAbHgwzhrNoSoXMeHBrM8tZn1uMQ+cO5SU+G78/fPtPHDuUGIjLDKvXnR6EvBFu5mR1o91O4qMdEr6wgxGeS6Q9okOpSBAnj6+rJgN/7zGaNsiojk9/VXKLWGN9qeBonL3+rRJceHkFZUzYVAPRg2IrhPUvd/LvHrR2UnAF23Kaivj8aWZzJmWSkp8N56ZOYJXvs5hycZ8rDYHFU73XPo9BxsGe6Vree2Dx5hi3WBsO/+6uWxLCByUY8PNWMxBxIYHk7m/jMSoUPaWVHLWCb2I84ziffPzEuBFVyIBX7Spx5dmekb0mbx5/RiKy52s2H6A4vJqzCZw1vi/eWrGli956vPnjfZjZ9zE66MvarI/e7mLsOBanrhkGFkFh4xFyr2jeLkoK7oyCfiiTc2Zlkp1zVYG94zE7nDy+NJMrDYHZpPCVaspLKs21pkFGFiUz1evzTaevzFxKDOv+gc1zViMxGyCIJO7NMLD/9nGFWP6ExN+OE0jxc5EV9eigK+Uegq4AHACVuB6rfVBpVQS8DOQ5dl1ndZ6tt+DiE7Lu9pUSbmT+at38HVWAfnF5ZhNipraw0UNSitdBNXW8I/Pn+fSrSuM7RNmv05+VEKT/USGmDlU5cJVC1eO6cu31iJOTenRILjLRVnR1bV0hL8MeEBr7VJK/QN4ALjP85hVaz2ihccXx7DHl2bWKXL2ywGH3/3G7drCQ8vnc7wtD4Dbpt/PZ0MnNKuPscmx3HfOUJ5d9gupvSO5bEx/+sSEMzU1gT4xYVLsTAgfLQr4WusvfZrrgEtbdjriWOHv5inviD4xKoS9JVVcltaXTbuKA646lVh6gAdXvsm07avJ796TWy56kC+OOyXgYiQAFgWxkSH0jAxhy55Sxg2MZeSAGKMsgqRthAisNXP4NwCLfNrJSqkfgFJgjtZ6dSv2JY6y+lUtl2Ts5qvtB1ifazf22by7mDI/FS1DXE7S13/AbeveR6F5dsJVzBtzCVXBIQ32DQJqOFwSwandc/QH9YzkgXOHNkjPSNpGiMCaDPhKqeVALz8P/Ulr/Ylnnz8BLuAdz2P7gP5a6yKl1CjgY6XUCVrr0voHUUqlA+kA/fv3/3U/hWhXdoeTcqeLu6YMZkZaP/5vWRYL1+0iJT68zn4Ngr3WnJW9jjlfvUb/kgL+O+RU/jb5RvZE9WzQh7f2TUxEMIWOaqOM8cj+0YQFB/HI9BNIiW+4WpWkbYQITGldvyL4ER5AqVnAbGCK1ro8wD5fA7/XWmc0dqy0tDSdkdHoLqID8KZN7poymHBLEAu/y/M7j95XSuFuHloxn9PyfiCrR38ennoL3w0Y3uw+xybHMm5gLLPGJ8uUSiHqUUpt1FqnNbVfS2fpnIP7Iu3pvsFeKRUP2LXWNUqpgcBgYEdL+hLtr365Yu/3M9L6Ue50sW5HIetzixuM7H1FVjm4a82/mbVpKeXBoTw09Rb+dfJ5TU6zNCn3AuKHqmqYMCiOR6afyPLMglb9+YToalqaw38RCAGWKfeFNu/0y9OAR5VSLtwp2Nlaa3vgw4iOpv7SgkCdi6H7SipZn1sMwO6ihh/slK7l0p9W8MdVC4grL+G94Wfz9GnXYA+Palb/tRoOVdUweUg8z8wcIYXNhGgFLZ2lMyjA9g+AD1pybHF0LVibx8osm1HkbNH3u5gwqAejk2KZt8rK5z/tM/Z11tZ97oi9WTy8/BVG7MtmY+JQrp/xMFt7+X2rAO7a9MHmIEorXfSLCTMWOgkLNvHMzBFS2EyIViJ32gq/KpzuqZSpvSNZ9P0u5q/OBaC6ppb1uXZiw4Nxf3g7rIejmPu+XsCMrcsp6BbL3dPu5ePUSY1Os4ywmDixTzTrc+2kxEcw/9o05q2y8uW2Av7vshFGvl4uxgrRchLwhaFuzt4dpCura/nSJ3eeaytz7+upQgkQXFPNtRuXcte37xLqcvLy2Et58ZSZOEIC5/a9osMt3HfOUF74KtsosPbkpcN5Uu7oEKLVScDvgvzdNAWH59Z/tf0Aew+60ypZBYfYX3p4Bs6BMmedY03M3cRDy+czyJ7PyoGjeHRKesBFw/3Zc7CSVb8cYNzAOGLCZfaNEG1JAn4X5O+mqdFJsazOtjGyf7Rx81SwCawHDvk9Rr+D+/nzV69xVvY6cmN6c/2lD7EyZXSzz6FfTBjnntgbFGzcWWwsVC5pGyHajgT8Lsh74XNqagJ3vvsDa3IKjYuliVGhxn7VtVDoqK7z3DBnJbeuW8It33+IyxTEE6dfxxtp03GagxvtMz7CwuTje7Imu5BeUaE8NWM4KfHdmLfKypqcIiYPiZcLskK0MQn4XVBshIWpqQnc+NYG8jxTKiND3W+FgtIAN1BpzbTtq3lw5RskHirko9RJPDHpOgoiezSrz4tH9mH2pEEsia+bSvKdfSM3VAnRtiTgd1GPL800gj1AaYWL4CCo9rMeyfEHdvDw8vmM3b2VrQkp3HnhH8joe0KTfQzrE8mBQ05MSnHOib39zqWX2TdCtB8J+F2M3eFkwdpc4iNDiAgx4aiqxQzkey7S+oqqOMS9q//FVZs/pyS0Gw+cfTuLhp1JbSN3ycaGm7GXu7CYAJSxGPkLX2XzzEx3tWxJ3QhxdEjA70L83T0L7qp3vky1NVzx4xf8/pu36V7l4O2Tz+O5CVdREhbZdB/l7qM5a2GnvZyR/aMIDgpizrRUGc0LcZRJwO+kvFMvvWu6Tk1N4KFPtrImp4iQIEVVjf+ieaN3b+Xh5fM54cAOvut/Eg9PvYWs+KQm++sXE8bkIfFkFZSxPteO2QQlFS7OPqG3BHkhOggJ+J2I7/z6BWtzmbsix6hRvzq70Jj66C/YJxwq5MGVbzL951XsiYx3rzo15NRG75Id1jeKULMJpRTrc+30iQnn7jOHGJ8iZOaNEB2LBPxjnG+Q970o6r1TFk8l+R93HyRIQf1Yb3FVc9OGj/jtd4sx19Ywd/zlvDL2UiosoTQlMsTMC1eOBKhzI5e32JnMvBGiY5GAf4zzBvnFGbv58/mpjE2O4Ytt++kXE0ZSXDgXjejDxp0HOVRVL1OvNVOs3/PnFa+RdHAf/zvuFB6ffCP50f7Wujmse6iZ0koXiVGhfGst4pWvc3jw/NQ6aRvJ1QvRMUnA70AClTxozNTUBF7+2orV5uCuRT8Y68du2nUQgAc/2kr9BM7Aonz+suJVJuVuJDuuH1fPfIw1ySc32VdYsIm/XnQie0sq3eUXSirJ3Of/TlwhRMcjAb8D+TU13z/dvJeDFe67Yf0tFu4b7LtVlXP72ve4IeNTKs0WHjvjJhaMnIYrqOm3QWSIezGSD3/Yw5vXj2FqagKPL81kzrTUZp2nEOLok4DfgTS35rvVVuYTbN0h3V9+3kvpWi7etpL7v36Lno5iFp10Jk+dfi2FETEB+zCbwFULPSNDOHCoisjQYFITo7jjjMHMW2VlRlo/3rx+zK/6OYUQR4cE/A6kubnvx5dmsjLLhtO1lRMSo7CYGi5C4nXSvmweXj6PUXu3s7n3cdx8yRx+TBzSjHNxB/oDh6qIjQhmb0klQ3pFsiHPLitPCXGMkoDfwdWfTz86KZZyZw0nJUaSV+TgW2uR3+fFOQ7y+28WctmWZRRFRPH78+7mgxPPQCtTo/2FBZswBykOHHLfITthUBy/O3OIUa/eW8JYplsKcexRWgfIAxwFaWlpOiMj42ifRocyb5WVv3++naS4cPKKyokKM/vN1XtZXNVctfkz7lnzb8KqK3lz1IW8cOrlHAqJaLKvfjFh9I8N51trEQNiw7no5ERmjU+WqZVCdHBKqY1a67Sm9pMRfgfkO6rfU+wO8t5CZ40F+99/s5Dbv1sMwDdJJ/PI1HSscQ1H4v7y/b26h/DWDe6cvPf6QEp8t1b6iYQQHYEE/A7EG+jLnS7mrshhccZurDZHk88bu+snFr37gNF+a+Q0Hp56S8C7ZE0cXo22R0QwV41LYtb4JGMkLxdjheicWhTwlVIPAzcD3mpcD2qtP/M89gBwI+7YcqfW+ouW9NUVeKdlpk8cSFJcOFabg24hQZRV1aCgwXz66IpSNj9/pdGuMIcw9rcLKA1tODI3c7hI2sCeEewsqqDSVUuYxVwn2AshOq/Gr+A1z3Na6xGef95gnwpcDpwAnAO8pJQKXFO3C7M7nMxbZcXucDIjrR93TRlExk67kcIpq3KPxesEe62Z++lTdYL9JVc/xfH3ftAg2FuC3OkaF+5UDsDYgT34710TSYmPYHdxBUsydrfhTyiE6CjaKqUzHXhPa10F5CqlcoAxwHdt1N8xxV/9m3U7iph1ShJvrc1rNE9//s+r+een/zDaT0+8mhfHXx5wf2cN7C+tIjrczEFP6eJQs4nlmQXMvzaN5ZkFMuNGiC6iNQL+7Uqpa4EM4F6tdTHQB1jns0++Z1sDSql0IB2gf//+rXA6HZvVVmYsLfjFtv0EB5kY1ieSlVk2vvmlkJoAs6b6HtzPmnk3Ge2sHv25YNbcJteS9bpwWB9iIiy4PysomUsvRBfUZMBXSi0H/FXU+hPwMvAY7ijyGPAMcAOHSzX68hvJtNbzgfngnpbZrLM+RtkdTtIXZhjpGm+9m6gw96/BX7A317hY8s59nLwvy9g2+eZ55Mb6/ftJ91ATpZXuu7BSe3Ujc38ZI/tHc/eZxxl5ervDSbglSEb2QnQxTQZ8rfXU5hxIKfUqsNTTzAd8o0lfYO8Rn10n4F5SMA/QVFTXYrU56B5qpm9MGFWuWvLt5VhM/i+l3Pj9R/x55etG+97z7uGDk6b43dc71TIpLpLuYWbW5BQRGeYO8BMHx9e5KCvVLIXomlo6S6e31nqfp3kxsNXz/afAv5VSzwKJwGDg+5b01dHVvyPWW/FyScZu5q7IBuDUlDgASitddapM2hzOOsc6YX8O/11wt9FeNmgM6ZfMafQuWe+8eldNLVrDteMGEBpsYtzAOGaNT2qln1IIcSxraQ7/SaXUCNzpmjzgFgCt9Tal1GIgE/dswN9qrWsCHqUT8L346l0z9pbTU5iamsBX2w8AmpsnDsRqKzMW9q4v3FnBN/Nuokd5ibEt7fa3Gy1y5uWtZnmwoprM/YfYX1qJ1eZg8pB4isuPvOyyEKLzaVHA11pf08hjfwX+2pLjH0u8+fCpqQkM67uXorIqnlv2CxXVNazPtQNQXfNLwGD/0PJ5XL/xP0b7mpmPsjp5ZKN9hppN9I8N45cDDgYnRLJp10GmHp/AjsIyBsZ3o1f3UM8fn8w6f4SEEF2T3GnbSmIjLMY0ywpnDfNX5wLQJ/rwUoGbd5U0eN7pOzayYMlDRvvNURfwyNRbGu0rMSqUfrFhrM8t5pcD7lH8nGmpRirJ+2njrimDOe24eKamJjBuoEy/FKKrk4DfSqy2MtIXZmC1OUiKCze27zlYaXzvW8E4vqyYDf88/AHJFhHN6emvUm4J83v8xCj3H45ardlbUsm0Yb0ZN7AHoI0CZymnu2+68q2r703heB8TQnRdEvBbwD0DJxdQZOTZsdocRIe7C515FxCpT+laXvvgMaZYNxjbzr9uLtsSAqdaIiwmpg1LZP7qHcY2b0kEf3fJyiwcIYQ/EvBbwD0DJwdwp276RocxakAMn/y412+wn7FlGU99PtdoP3bGTbw++qIm+/nNyH7MnpRCmMVEhbOWMEuQEezlBiohRHNJwP8VNu0s5g/v/8ifz0/l2nH9WbQh30jd7DlY0WD/gUX5fPXabKO9MXEoM6/6BzWm5pUXCg02ERth4Z4z665U1dwlEYUQAiTgHzG7w8l1b35PaaWLP76/BZeuparm8HDe915Zi6ua/751J4OLDqddJsx+nfyohIDHN+HO9QeZFDW17qOFWfz/miR1I4Q4EhLwm8G32NmCtXmUVrqLkBU6qqgNUAzirjX/5p5v/220f3vhffz3+IlN9lULJMWF8/p1o/l08x5AyY1TQohWIQG/GRaszWXuihz2FFew/OcCY7u/YD8y/2c+fOcPRvuj1EncM+3egIuRhJnBWxzTu3zhmORYUuK7MWt8spQuFkK0Ggn4fviO6AFWZxcC8L9t+43FvevrXlnGhhevIaSmGgCXMjHqjncoCYtstK8KF/TuHoo5SBETbmbLnkPsKXZfB5CLskKI1iQBvx67w8m9izezMsvGu9/vIjbCYlS1LPQX7LXm6c/+j0u3rjA2zbzyCb7vd6Lf40dYTJw/LJG8Qgff5xUbi5MDJEbHGjdRgVyUFUK0Lgn49SxYm8fKLBshQYq8onIjGCvq3jgFcPYva5n30d+M9vOnXMazpwWsNgGAw1mL7VAVf//NMJZnFjA1NYFF3+8mc18Jj0w/sc7C4XJRVgjRmiTgcziFMzoplg835QNQVVM3Qe/b6l1q47uXrzfaO2ISOfeGF6kyN12YLCrMzMosG+MGFhjB/MHzj2/5DyGEEE3osgG//sybuSuy6R5qNmbg+BNUW8O/332QsfnbjG1Tb3yJnB6Nr9RlUu4LvElx4Tw7cwQb8uxMTU1g3iqrVLAUQrSbLhvw/2/ZLyxct5OF3+VxakoPACPYKxouz3Xtxv/w6PJ5Rvv+s2/nvRHnNDhu/ecqBWelJvC/bQVMH5HIyAExjBwQw7xVVrkgK4RoV1024H+T7S4XvOdgJR9t3lPnMd+APfRALv978w6jvSp5JNfNeDjgYiThFhMOp8+NWBq27S1lwqAeVDhrsDucRmVNkAuyQoj20yUD/qadxVT73B1bXdNwQn1odSUr599C77IiY9uY2xZwIDKu0WOn9IxkS34JQSaoqYXYiGAmD+nJwnU7WZNTSFy3EG45PUUuyAoh2l2nDfi+a8l6ywd7q1su/G4nxeXVAZ97/8o3mP39h0b7+ksfYmXK6Cb7nDCoB49MP4GHPtnGmpxCxibHMm5gHBeOSCQmwgJoGdELIY6aThvwvXfHurkz6xt3HmRNTqHPlrrG523m34vmGO23Tz6PP591W7P6G5scy/NXnExshIVRA6JZk1NIcJBi7opswi1B3HPmcS3+mYQQoiU6XcD3zr6pqHanbMYmx7ByewFb9pQCEGEJwuGsqRPsY8tL2PTCVUa7NCSC8be+SVlIOM3hHsnHGu1Z45MJt5jrLGguhBBHW6cL+N5yBElx4Vx7ygCy9h8ygj2Aw+mzlrrWvPTx3znvl7XGpguvfZYtvZsejQeboLrWPdVyeN8o5q7IIdxibpCfl5WmhBAdRYsCvlJqEeAt0h4NHNRaj1BKJQE/A1mex9ZprWc3PELr8b15KiY8mLyicgo35VNWVeN3/+nbVjJ36TNG+4nTr+OVcZc2q6/I0CDOGNKTH/NLyCsqJ8xi5oFzh8pIXgjRobUo4GutL/N+r5R6BvBdpduqtR7RkuMfCe/Ivlf3EOOCrL9gP6B4L6vmpxvtLb0GccnVT+MKCvxSBAE1QHyEhYtH9SUs2GRcH5g8JJ5Z45Pk5ikhRIfXKikdpZQCZgJntMbxfo0Zaf1Yt6OIlVk2v48H11Tz0du/58QCq7HttPRX2RXTu9Hj+hY3u3LcAO458zjsDqfnUSXBXghxzGitHP5EoEBrne2zLVkp9QNQCszRWq9upb78io2wMGdaKut2rDYu2Hrdum4J961aYLTvmnYvn5wwucljdg818/AFJ/Damh2k9o4yFiLxt9ygEEJ0dE0GfKXUcqCXn4f+pLX+xPP9FcC7Po/tA/prrYuUUqOAj5VSJ2itS+sfRCmVDqQD9O/feE2axlhtZaQvzKgT7IfvzeKTt+812kuHTOD26fcFXIzEBKTEh7PTXoGzRlNa6eLlVVbW59oZNSDGqL0jI3ohxLGoyYCvtZ7a2ONKKTNwCTDK5zlVQJXn+41KKStwHJDh5/jzgfkAaWlpARYMbNrd727CanMA0K2qnO9emkWk8/CC4iPveAd7eFSjx+jZPYRsmzt94y2kVl5VzeQh8VQ4a5m7QmrfCCGOXf4LwhyZqcB2rXW+d4NSKl4pFeT5fiAwGNjRCn0FlFdUDlrz1y9eZOv/zTSC/RWX/5Wk+5YGDPY9Itx/86LCzOwvdS9wMjY5lotG9AEgPCSYlVk2wixBMhNHCHFMa40c/uXUTecAnAY8qpRy4Z7gMltrbW+FvgI6O3cjT799+C7ZV8b+hicmXe93X7MJXJ7Mz3nD+pC1v5T1ucUkRoWyt6SScQNjmTU+mT4xYXVunpJUjhDiWKa0/tVZlFaXlpamMzIaZH2advAgxMQAsCcynik3v0xlcGiTT+sXE8ZbN4zhoU+2sianiH4xYVwysq/MvBFCHFOUUhu11mlN7dc57rSNiuLRi37H2phktvdM9ruLAgb3jGBcSg9jRL+7uILlmQU8Mv1E0hdmYLU5CLcESbAXQnRKnWOED3y9/QDXv7WhQUE0X9eOG0BMRLB7Jo+GMIupTiVNmYUjhDgWda0RPrDWWlQn2AcHKaprNCYOLz6+MusAu4vdF3MnD4nnmZkjjOAu9emFEJ1da8zS6RAydta9Jlxdo4mNCOb9W8czNtldyXLykJ7cNWUQEwb1YGWWjSUZuwF3HZ55q6w+d9AKIUTn02lG+MFBdW+migozY3dUsyHPzstXj6qTrvFN38DhOjwgc+yFEJ1Xp8nhb9pZzD2LNhMVZmbLnlLSJw4krpulWTl5yd8LIY5lzc3hd5qUzoY8Ozvt5UwemsAD5w7lsjHNv0HKm7+XYC+E6Mw6TUpnamoC63YUceGIRFLiuzFvlVXSNEII4aPTBPzlmQWszLIxbmABKad3M/LzUgpBCCHcOk3Arx/gZZqlEELU1Wly+EIIIRrXaQK+d2qld269EEKIujptSkcIIURdnSbgS85eCCEa12lSOkIIIRonAV8IIboICfhCCNFFSMAXQoguQgK+EEJ0ERLwhRCii5CAL4QQXUSHqoevlLIBO1twiB5AYSudTmuS8zoycl5HRs7ryHTG8xqgtY5vaqcOFfBbSimV0ZxFANqbnNeRkfM6MnJeR6Yrn5ekdIQQoouQgC+EEF1EZwv484/2CQQg53Vk5LyOjJzXkemy59WpcvhCCCEC62wjfCGEEAEcUwFfKTVDKbVNKVWrlEqr99gDSqkcpVSWUursAM9PVkqtV0plK6UWKaUsbXSei5RSmz3/8pRSmwPsl6eU+smzX0ZbnEu9/h5WSu3xObfzAux3jud1zFFK3d8O5/WUUmq7UmqLUuojpVR0gP3a5fVq6udXSoV4fsc5nvdTUludi0+f/ZRSK5VSP3v+D9zlZ59JSqkSn9/vX9r6vDz9Nvp7UW7Pe16vLUqpke1wTkN8XofNSqlSpdTd9fZpl9dLKfWGUuqAUmqrz7ZYpdQyTxWeT4IAAAT8SURBVCxappSKCfDcWZ59spVSs1p8MlrrY+YfcDwwBPgaSPPZngr8CIQAyYAVCPLz/MXA5Z7vXwFubYdzfgb4S4DH8oAe7fj6PQz8vol9gjyv30DA4nldU9v4vM4CzJ7v/wH842i9Xs35+YHbgFc8318OLGqH311vYKTn+0jgFz/nNQlY2l7vp+b+XoDzgM8BBYwD1rfz+QUB+3HPVW/31ws4DRgJbPXZ9iRwv+f7+/2954FYYIfna4zn+5iWnMsxNcLXWv+stc7y89B04D2tdZXWOhfIAcb47qCUUsAZwPueTQuAi9ryfD19zgTebct+WtkYIEdrvUNr7QTew/36thmt9Zdaa5enuQ7o25b9NaE5P/903O8fcL+fpnh+121Ga71Pa73J8/0h4GegT1v22YqmAwu12zogWinVux37nwJYtdYtuanzV9NafwPY6232fQ8FikVnA8u01natdTGwDDinJedyTAX8RvQBfBezzafhf4Y44KBPYPG3T2ubCBRorbMDPK6BL5VSG5VS6W18Ll63ez5WvxHgY2RzXsu2dAPu0aA/7fF6NefnN/bxvJ9KcL+/2oUnhXQysN7Pw6copX5USn2ulDqhnU6pqd/L0X5PXU7gQdfReL0AErTW+8D9xxzo6WefVn/dOtwSh0qp5UAvPw/9SWv9SaCn+dlWf/pRc/Zptmae5xU0Pro/VWu9VynVE1imlNruGQ38ao2dF/Ay8Bjun/sx3OmmG+ofws9zWzyVqzmvl1LqT4ALeCfAYVr99fJ3qn62tel76UgopboBHwB3a61L6z28CXfaosxzfeZjYHA7nFZTv5ej+XpZgAuBB/w8fLRer+Zq9detwwV8rfXUX/G0fMB39fK+wN56+xTi/ihp9ozK/O3TbE2dp1LKDFwCjGrkGHs9Xw8opT7CnU5oUQBr7uunlHoVWOrnoea8lq1+Xp4LUtOAKdqTwPRzjFZ/vfxozs/v3Sff83uOouFH9lanlArGHezf0Vp/WP9x3z8AWuvPlFIvKaV6aK3btG5MM34vbfKeaqZzgU1a64L6Dxyt18ujQCnVW2u9z5PeOuBnn3zc1xm8+uK+fvmrdZaUzqfA5Z7ZE8m4/0p/77uDJ4isBC71bJoFBPrE0BqmAtu11vn+HlRKRSilIr3f475wudXfvq2lXt704gD9bQAGK/eMJgvuj8OftvF5nQPcB1yotS4PsE97vV7N+fk/xf3+Aff76atAf6Rai+cawevAz1rrZwPs08t7LUEpNQb3/++iNj6v5vxePgWu9czWGQeUeNMZ7SDgp+yj8Xr58H0PBYpFXwBnKaViPOnXszzbfr22vkLdmv9wB6l8oAooAL7weexPuGdXZAHn+mz/DEj0fD8Q9x+CHGAJENKG5/oWMLvetkTgM59z+dHzbxvu1EZbv35vAz8BWzxvuN71z8vTPg/3LBBrO51XDu5c5WbPv1fqn1d7vl7+fn7gUdx/kABCPe+fHM/7aWA7vEYTcH+c3+LzOp0HzPa+z4DbPa/Nj7gvfo9vh/Py+3upd14K+Kfn9fwJnxl2bXxu4bgDeJTPtnZ/vXD/wdkHVHvi1424r/msALI9X2M9+6bB/7dvhzYAgDAUBbsqzMVoDINB4BA4/p1vgmieaEKNY7bvPZtV1V7f4qctQIhfTjoAXAg+QAjBBwgh+AAhBB8ghOADhBB8gBCCDxBiAXdbGCReXiNRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Optional\n",
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(X,y, s=1)\n",
    "plt.plot(X,true_w * X + true_b, 'red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T03:38:01.470849Z",
     "start_time": "2019-07-13T03:38:01.467199Z"
    }
   },
   "source": [
    "# 3. Gradient Descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T05:54:07.874858Z",
     "start_time": "2019-07-13T05:54:07.870031Z"
    }
   },
   "outputs": [],
   "source": [
    "def gradient_descent(X, y, w, b, learning_rate):\n",
    "    dw = -2 * np.sum(X * (y - w * X - b)) # ∂e/∂w\n",
    "    db = -2 * np.sum(y - w * X - b)       # ∂e/∂b\n",
    "    w_new = w - learning_rate * dw        # minus sign since we are minizing e\n",
    "    b_new = b - learning_rate * db\n",
    "    return w_new, b_new\n",
    "\n",
    "def get_loss(X,y,w,b):\n",
    "    return (y - w * X - b).T @ (y - w * X - b)   # square loss, \n",
    "    # .T and @ denote transpose and matrix multiplication resp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T06:00:56.923946Z",
     "start_time": "2019-07-13T06:00:56.867108Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "103801\n",
      "91779\n",
      "90183\n",
      "89971\n",
      "89943\n",
      "89939\n",
      "89939\n",
      "89939\n",
      "89939\n",
      "89939\n",
      "y = 7.59 x - 3.26\n"
     ]
    }
   ],
   "source": [
    "learning_rate = 0.000001\n",
    "max_epoch = 500\n",
    "\n",
    "w, b = -1,0\n",
    "\n",
    "for epoch in range(1,max_epoch+1):\n",
    "    w,b = gradient_descent(X, y, w, b, learning_rate)\n",
    "    \n",
    "    if epoch % 50 == 0:\n",
    "        print(f'{get_loss(X,y,w,b):.0f}')\n",
    "\n",
    "if b > 0:\n",
    "    print(f'y = {w:.2f} x + {b:.2f}')\n",
    "else:\n",
    "    print(f'y = {w:.2f} x - {-b:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Verify using sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-07-13T06:00:37.001434Z",
     "start_time": "2019-07-13T06:00:36.993307Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R^2 score: 0.995308982100681\n",
      "y = 7.59 x - 3.26\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "reg = LinearRegression().fit(X.reshape(-1, 1), y.reshape(-1, 1))\n",
    "print(f'R^2 score: {reg.score(X.reshape(-1, 1), y.reshape(-1, 1))}')\n",
    "\n",
    "w, b = [np.asscalar(v) for v in [reg.coef_, reg.intercept_]]\n",
    "if b > 0:\n",
    "    print(f'y = {w:.2f} x + {b:.2f}')\n",
    "else:\n",
    "    print(f'y = {w:.2f} x - {-b:.2f}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
